#### 0.准备好环境 ####
# Paper02-Sample_02_paper_成人非酒精性脂肪性肝病患者血清维生素 C 水平与肝纤维化的关系
# 本文件是测验用的
# 引入包
library(haven)
library(plyr)
library(dplyr)
#加权分析用
library(survey)
#tableby
library(arsenal)
library(gtsummary)

# 4.3节用到包
library(tidyverse)
library(gtsummary)
library(tidyr) # drop_na 函数，用于快速去掉 NA

setwd("G:/BaiduNetdiskDownload/NHANES/")
#### 1. 定位数据模块和变量，获取源数据 ####
##### 1.1 DEMO-人口学数据提取 #####
demo.j <- read_xpt("2017-2018/Demographics/demo_j.xpt")
# colnames(demo.j)

### 2 提取研究所需要的变量
#  提取这些
# 年龄 RIDAGEYR
# 性别 RIAGENDR
# 种族 RIDRETH3
# 教育水平 DMDEDUC2
demo.data <- demo.j[,c('SEQN', 'RIDAGEYR', 'RIAGENDR', 'RIDRETH3', 'DMDEDUC2')]

##### 1.2 PAQ-运动数据提取 ##########
### 1.提取 Component文件
# https://wwwn.cdc.gov/Nchs/Nhanes/2017-2018/PAQ_J.htm
paq.j <- read_xpt("2017-2018/Questionnaire/paq_j.xpt")
### 2.提取研究所需要的变量
# 2023年2月22日14:54:40 这里不需要转换 直接拿取就可以了
# paq.data.file <- dplyr::bind_rows(list(paq.j))

paq.data <- paq.j[,c('SEQN', 'PAQ650','PAQ665')]


##### 1.3 BMI 数据提取 ##########
# https://wwwn.cdc.gov/Nchs/Nhanes/2017-2018/BMX_J.htm
### 1.提取 Component文件
bmx.j <- read_xpt("2017-2018/Examination/bmx_j.xpt") #注意是 Examination 的类别

### 2.提取研究所需要的变量 BMI-BMXBMI; 腰围-BMXWAIST
bmx.data <- bmx.j[,c('SEQN', 'BMXBMI')]

##### 1.4 SMQ-吸烟数据提取 ######
# https://wwwn.cdc.gov/Nchs/Nhanes/2017-2018/SMQ_J.htm
### 1.1 提取 Component文件
smq.j <- read_xpt("2017-2018/Questionnaire/smq_j.xpt")

### 2. 提取研究所需要的变量
# 是否吸烟至少100支-SMQ020; 现在是否吸烟-SMQ040; 多久前开始戒烟-SMQ050Q;
# 多久前开始戒烟的时间单位（天、周、月、年）-SMQ050U;
#smq.data.file <- dplyr::bind_rows(list(smq.g, smq.h))
# colnames(smq.j)
smq.data <- smq.j[,c('SEQN',"SMQ020", 'SMQ040', 'SMQ050Q')]

##### 1.5 糖尿病诊断  #####
# Participants with diabetes were identified as having any of the following: 
# (a) hemoglobin A1C concentration >= 6.5% ***or*** a fasting plasma glucose level >= 126 mg/dL [21]; 
# (b) for those who responded ‘yes’ to the question: ‘Doctor told you have diabetes?’ ***or*** ‘Taking insulin now?’.

ghb.j <- read_xpt("2017-2018/Laboratory/ghb_j.xpt")
ghb.data <- ghb.j[,c('SEQN','LBXGH')]

##### 1.5.2 空腹血糖（mg/dl）##### 
# https://wwwn.cdc.gov/Nchs/Nhanes/2017-2018/GLU_J.htm
glu.j <- read_xpt("2017-2018/Laboratory/glu_j.xpt")
# colnames(glu.j)
glu.data <- glu.j[,c('SEQN','LBXGLU')]

##### 1.5.3 是否有医生告知您患有糖尿病、现在是否使用胰岛素 ##### 
diq.j <- read_xpt("2017-2018/Questionnaire/diq_j.xpt")
# colnames(diq.j)
diq.data <- diq.j[,c('SEQN','DIQ010','DIQ050')]

#合并提取的四个维度的数据
diabetes <- plyr::join_all(list(ghb.data,glu.data,diq.data),by = "SEQN")
diabetes.a.index <- ifelse(diabetes$LBXGH>=6.5 | diabetes$LBXGLU >=126,1,NA)
diabetes.b.index <- ifelse(diabetes$DIQ010==1 | diabetes$DIQ050 ==1,1,NA)
# length(diabetes$LBXGH)
# length(diabetes$LBXGLU)
#这两个指标合并一起
diabetes.index <- ifelse(diabetes.a.index==1 | diabetes.b.index==1,1,0)
# 附件一起
diabetes$diabetes.index <- diabetes.index

length(diabetes)



##### 1.6 丙氨酸转氨酶  ##########
# https://wwwn.cdc.gov/Nchs/Nhanes/2017-2018/BIOPRO_J.htm
alt.j <- read_xpt("2017-2018/Laboratory/biopro_j.xpt")
# colnames(alt.j)
alt.data <- alt.j[,c('SEQN',"LBXSATSI")]

##### 1.7 血清维生素C的含量 #####
# https://wwwn.cdc.gov/Nchs/Nhanes/2017-2018/VIC_J.htm
vic.j <- read_xpt("2017-2018/Laboratory/vic_j.xpt")
# colnames(vic.j)
vic.data <- vic.j[,c('SEQN',"LBXVIC","LBDVICSI")]

##### 1.8 高密度脂蛋白 #####
# https://wwwn.cdc.gov/Nchs/Nhanes/2017-2018/HDL_J.htm
hdl.j <- read_xpt("2017-2018/Laboratory/hdl_j.xpt")
# colnames(hdl.j)
hdl.data <- hdl.j[,c('SEQN',"LBDHDDSI","LBDHDD")]


##### 1.9 弹性成像检查  #####
# https://wwwn.cdc.gov/Nchs/Nhanes/2017-2018/LUX_J.htm
# 弹性成像检查状态-LUAXSTAT
# 可控衰减参数中位数-LUXCAPM
# 中位硬度-LUXSMED
lux.j <- read_xpt("2017-2018/Examination/lux_j.xpt")
# colnames(lux.j)
lux.data <- lux.j[,c('SEQN',"LUXCAPM","LUXSMED","LUAXSTAT")]

##### 1.10 饮酒数据提取 #####
# https://wwwn.cdc.gov/Nchs/Nhanes/2017-2018/DR1TOT_J.htm
dr1tot.j <- read_xpt('2017-2018/Dietary/dr1tot_j.xpt')
# colnames(dr1tot.j)
dr1tot.data <- dr1tot.j[,c('SEQN',"DR1TALCO")]

dr2tot.j <- read_xpt('2017-2018/Dietary/dr2tot_j.xpt')
dr2tot.j$DR2TALCO
# colnames(dr2tot.j)
dr2tot.data <- dr2tot.j[,c('SEQN',"DR2TALCO")]


alco.data <- merge(dr1tot.data, dr2tot.data)
# dim(alco.data)
# #View(alco.data)

##### 1.11 其他原因导致的肝脏疾病-viral hepatitis infection #####
# 乙型肝炎表面抗原
viral.hepbd.j <- read_xpt('2017-2018/Laboratory/hepbd_j.xpt')
viral.hepbd.data <- viral.hepbd.j[,c('SEQN',"LBDHBG")]

# 丙肝抗体\丙肝RNA
viral.hepc.j <- read_xpt('2017-2018/Laboratory/hepc_j.xpt')
viral.hepc.data <- viral.hepc.j[,c('SEQN',"LBXHCR", "LBDHCI")]

# 肝癌、自身免疫性肝炎
mcq.j <- read_xpt('2017-2018/Questionnaire/mcq_j.xpt')
liver.disease.data <- mcq.j[,c('SEQN',"MCQ510E", "MCQ230A", "MCQ230B", "MCQ230C")]


#### 2 提取分析相关变量（权重等，暂时为复现paper结果而提取） ####
##### 2.1 权重变量 ##### 
# 找到权重变量这块我是看不懂
# DEMO->MEC-Dietary, GLU(空腹血糖)
# 饮食中的权重：https://wwwn.cdc.gov/Nchs/Nhanes/2017-2018/DR1TOT_J.htm
# WTDRD1, WTDRD2

# GLU: https://wwwn.cdc.gov/Nchs/Nhanes/2017-2018/GLU_J.htm
# WTSAF2YR

# 权重的计算方式，需要按 Cycle 进行平均，这个仅有1个 Cycle，无需做处理
weight.data <- demo.j[,c('SEQN', 'WTMEC2YR')] 

##### 2.2 复杂抽样的其他变量-DEMO #####
survey.design.data <- demo.j[,c('SEQN', 'SDMVPSU', 'SDMVSTRA')]

#### 3.合并上述所有数据（把列拼接起来）-Output #### 
# 【4个# 用作第一层级】

##### 3.1 合并步骤1 & 2中提取的数据 #####
output <- plyr::join_all(list(demo.data, paq.data, bmx.data, smq.data,
                              ghb.data, glu.data, diq.data, diabetes,
                              alt.data, vic.data, hdl.data, lux.data, alco.data,
                              viral.hepbd.data, viral.hepc.data,liver.disease.data,
                              weight.data, survey.design.data), by='SEQN', type='full')

# #View(output) # 跑通到这一步，得到这篇文章所用到的全部原始数据
#### 4.根据文章的筛选策略进行筛选 ####
# 原始筛选策略：
#####  4.1缺少血清维生素 C 数据(n = 2514) #####
dim(output)

vc.exclude.index <- which(is.na(output$LBDVICSI))
length(vc.exclude.index)

#删除这些没有数据的 有- 是排除这些人没有- 是提取这些人
data.vc.exist <- output[-vc.exclude.index, ]
#dim(data.vc.exist) #6740

##### 4.2 排除 LUX 数据缺失以及不是 NAFLD 的人员 #####
# 弹性成像检查状态-LUAXSTAT
# 可控衰减参数中位数-LUXCAPM
# 中位硬度-LUXSMED
table(output$LUAXSTAT)
# complete ineligible   not done    partial 
# 5494        258        156        493 
# 这块排除需要仔细阅读原文并且有点难啊这块什么意思
cap.exclude.index <- which(data.vc.exist$LUAXSTAT != 1 ## 弹性成像检查状态-LUAXSTAT
                           | is.na(data.vc.exist$LUXCAPM)  #可控衰减参数中位数-LUXCAPM
                           | is.na(data.vc.exist$LUXSMED)  # 中位硬度-LUXSMED
                           | data.vc.exist$LUXCAPM <248) # 可控衰减参数中位数-LUXCAPM
length(cap.exclude.index) # 3945

data.vc.exist.cap.exist <- data.vc.exist[-cap.exclude.index,]
dim(data.vc.exist.cap.exist) #[1] 2795

##### 4.3 排除其他原因导致的肝脏疾病 & 重要协变量的缺失 #####
# 乙型肝炎表面抗原-LBDHBG, # 丙肝抗体\丙肝RNA-LBDHCI, LBXHCR, # 肝癌MCQ230A-C、自身免疫性肝炎-MCQ510E
# 肝癌 MCQ230A-C
liver.cancer.index <- which(data.vc.exist.cap.exist$MCQ230A == 22|data.vc.exist.cap.exist$MCQ230B == 22|
                              data.vc.exist.cap.exist$MCQ230C == 22)

# 自身免疫性肝炎-Autoimmune hepatitis
autoimm.hepa.index <- which(data.vc.exist.cap.exist$MCQ510E == 5)

# 丙肝抗体\丙肝RNA-LBDHCI, LBXHCR
hepc.index <- which(data.vc.exist.cap.exist$LBDHCI == 1|data.vc.exist.cap.exist$LBXHCR == 1)

# 乙型肝炎表面抗原-LBDHBG
hepbd.index <- which(data.vc.exist.cap.exist$LBDHBG == 1)

# 过度饮酒
# 先计算2天的平均值，其中第二天为 NA 的，平均值也是 NA
total.alco <- apply(data.vc.exist.cap.exist[,c('DR1TALCO', 'DR2TALCO')], 1, mean)
data.vc.exist.cap.exist$total.alco <- total.alco
#View(data.vc.exist.cap.exist[,c('DR1TALCO', 'DR2TALCO', 'total.alco')])

# 把第二天为 NA 的值计算出来，用第一天的值作为平均值
day.2.na.index <- which(is.na(data.vc.exist.cap.exist$DR2TALCO))
# 看下替换前的结果
 #View(data.vc.exist.cap.exist[day.2.na.index,c('DR1TALCO', 'DR2TALCO', 'total.alco')])


data.vc.exist.cap.exist$total.alco[day.2.na.index] <- data.vc.exist.cap.exist$DR1TALCO[day.2.na.index]
# 看下替换后的结果
#View(data.vc.exist.cap.exist[day.2.na.index,c('DR1TALCO', 'DR2TALCO', 'total.alco')])

# 过度饮酒的定义：Excessive alcohol consumption 
excessive.alco.male <- ifelse(data.vc.exist.cap.exist$RIAGENDR == 1 & data.vc.exist.cap.exist$total.alco > 20, 1, 0)
excessive.alco.female <- ifelse(data.vc.exist.cap.exist$RIAGENDR == 2 & data.vc.exist.cap.exist$total.alco > 10, 1, 0)
table(excessive.alco.female)

excessive.alco.index <- which(excessive.alco.male == 1 | excessive.alco.female == 1, 1, 0)

# 汇总上述指标
other.cause.index <- c(autoimm.hepa.index, 
                       hepc.index, hepbd.index,
                       liver.cancer.index, 
                       excessive.alco.index)
# 最后结果的汇总
data.vc.exist.cap.exist.non.other.cause <- data.vc.exist.cap.exist[-other.cause.index,]
dim(data.vc.exist.cap.exist.non.other.cause) #2374  

# covariate 上的缺失数据 也就是协变量缺失的人
paper.data <- subset.data.frame(data.vc.exist.cap.exist.non.other.cause,
                                (!is.na(LBDVICSI)) & # 维C，0，之前排除过缺失
                                  (!is.na(LBXSATSI)) & # 丙氨酸转氨酶, 28
                                  (!is.na(LBDHDDSI)) & # 高密度脂蛋白, 15
                                  (!is.na(LUXCAPM)) & # 可控衰减参数中位数-LUXCAPM，0，之前排除过
                                  (!is.na(LBXGH)) & # 糖化血红蛋白A1C(%), 2
                                  (!is.na(total.alco)) & # 24小时回忆的饮酒量
                                  (!is.na(RIDAGEYR)) & # 年龄
                                  (!is.na(RIAGENDR)) & # 性别
                                  (!is.na(RIDRETH3)) & #种族
                                  (!is.na(DMDEDUC2))  # 教育程度
)


dim(paper.data) # [1] 1926   41

#A median liver stiffness of 8.2 kPa was used to define cases of significant fibrosis (F2) in this study 
# 用于区分定义严重纤维化 分值



#吸烟原文判断 如果
# Participants were asked whether they had ever smoked 100 cigarettes in their lifetime 
# and whether they smoked currently to identify current
# and former smokers. Participants were defined as former
# smokers if they did not smoke currently but had ever
# smoked 100 cigarettes in the past
# 受访者问到一生中是不是吸过100 只烟 以及当前是不是吸烟来判断现在吸烟和以前吸烟 下面是解释
# 被定义为曾经吸烟者的情况是现在不吸烟但是在过去曾经吸过100只烟

# 一点也不吸烟 没说怎么判断 用我的方法得出的数量与论文不一致那么就是去除下面两种情况的都是从不吸烟
dim(paper.data[which(paper.data$SMQ040==3),])
# 当前吸烟
dim(paper.data[which(paper.data$SMQ040==1|paper.data$SMQ040==2),])
#曾经吸烟
dim(paper.data[which(paper.data$SMQ040==3 ),])


 paper.data$smoke_group = ifelse((paper.data$SMQ040==1|paper.data$SMQ040==2) , "Current",
                          ifelse((paper.data$SMQ040==3) , "Former" ,'Never'))
 paper.data$smoke_group[which(is.na(paper.data$smoke_group))] <-'Never'
 
 dim(paper.data[which(paper.data$smoke_group=='Never' ),])
 # nasmq040<-which(is.na(paper.data$SMQ040))
 # naSMQ050Q<-which(is.na(paper.data$SMQ050Q))
 # paper.data$smoke_group[nasmq040] <-"Never smoker"
 # paper.data$smoke_group[naSMQ050Q] <-"Never smoker"










